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Abstract 



This paper develops a new multivariate control charting method for vector autocor- 
related and serially correlated processes. The main idea is to propose a Bayesian mul- 
tivariate local level model, which is a generalization of the Shewhart-Dcming model for 
autocorrelated processes, in order to provide the predictive error distribution of the pro- 
cess and then to apply a univariate modified EWMA control chart to the logarithm of the 

■ Bayes' factors of the predictive error density versus the target error density. The resulting 
chart is proposed as capable to deal with both the non-normality and the autocorrela- 
tion structure of the log Bayes' factors. The new control charting scheme is general in 
application and it has the advantage to control simultaneously not only the process mean 
vector and the dispersion covariance matrix, but also the entire target distribution of the 

[ process. Two examples of London metal exchange data and of production time series data 

■ illustrate the capabilities of the new control chart. 
Some key words: time series, SPC, multivariate control chart, state space model, 



EWMA. 

OO ' 



1 Introduction 



In the last decades multivariate Statistical Process Control (SPC) has received considerable 
attention, since in practice many processes are observed in a vector form (Montgomery 1 ). 
Univariate control charts have been extensively discussed in the literature (Montgomery 1 , 
Box and Lucehno 2 , del Castilo 3 ) and many efforts have been devoted to upgrading the con- 
trol charts for: (a) cases of correlated univariate processes; and (b) cases of multivariate 
uncorrelated processes. 

Multivariate control charting has been discussed in many studies, e.g. Tracy et al. , Liu 5 , 
Kourti and MacGregor 6 , Mason et al. 7 , Vargas 8 , Ye et al. 9 and Pan 10 among many others. 
Review papers on multivariate control charts include Lowry and Montgomery 11 , Sullivan 
and Woodall 12 , Montgomery and Woodall 13 , Bersimis et al. 14 and Yeh et al. 15 . Most of 
the current research has been focused on the Hotelling's T 2 control chart and the multivari- 
ate EWMA control chart for controlling the process mean. Yeh et al. 1 ®, Surtihadi et al. 17 , 
Cheng and Thaga 18 and Costa and Rahim 19 propose and study multivariate EWMA and 
CUSUM control charts to control the dispersion of a multivariate process. As stated before 
univariate control charts for autocorrelated processes have been discussed in the literature 
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(Montgomery 1 , Box and Lucehno 2 ), however, for multivariate processes the general focus has 
been placed to uncorrelated processes. Dyer et al. 20 , Jiang 21 , Kalgonda and Kulkarni 22 and 
Noorossana and Vaghefi 23 consider multivariate control charting for autocorrelated processes 
based on autoregressive-moving-average (ARMA) time series models and the T 2 and mul- 
tivariate CUSUM control charts are illustrated. Pan and Jarrett 24 build a multivariate T 2 
control chart for the forecast errors of the process. They consider a state-space approach for 
modelling the underlying process and they point out that the problem of monitoring multi- 
variate processes is a problem of multivariate time series forecasting as well as a problem of 
control charting. Some forms of Bayesian control charts, known also as adaptive or dynamic 
control charts, are discussed in Tagaras 25 , Tagaras and Nikolaidis 26 , de Magalhaes et al. 27 
and in references therein. Adaptive control charts offer the flexibility and versatility to dy- 
namically change the sampling size and the sampling interval of a Shewhart control chart, 
but they are disadvantaged in that the complexity is increased and usually the modeller has 
to resort to Monte Carlo simulation. 

Our aim in this paper is to construct a multivariate control chart for autocorrelated 
processes in such a way that the scheme will be capable to monitor the process mean vector 
only, the process dispersion covariance matrix only, or both the process mean vector and 
the process dispersion covariance matrix. We propose a new control chart based on the 
theory of sequential Bayes' factors (West and Harrison 28 ). First we fit a local level model 
to the multivariate process and then we apply a univariate modified EWMA control chart 
to the logarithm of the Bayes' factor to monitor the dispersion of the predictive distribution 
of the data from the target distribution. Our model makes use of a generalization of the 
Shewhart-Deming model for multivariate autocorrelated processes (Deming 29 , del Castilo 3 , 
Triantafyllopoulos et a/. 30 ). 

Section [2] gives the necessary time series background. The proposed control chart is dis- 
cussed in detail in Section [3l In Sections 2] and [5] two examples, consisting of data from the 
London metal exchange and from a production of a plastic mould, illustrate the methodol- 
ogy and give light to the design and implementation of the new control chart. Concluding 
comments are given in Section [6] and the appendix details a proof of an argument in Section 

13 

2 Background 

The conventional control charts are based on the Shewhart-Deming model, e.g. for a p x 1 
process vector y t this model sets 

y t =fi + e t , e t ~A/" p (0,E), (1) 

where [i is the process mean vector and £ is the process dispersion covariance matrix, known 
also as the measurement covariance matrix. Here A/" p (0, S) indicates the p-dimensional normal 
distribution with mean vector zero and covariance matrix S. The measurement drift sequence 
{et} is assumed uncorrelated and this makes the generating process {yt} an uncorrelated 
sequence too. In this paper we extend the above model by considering equation ([I]), but 
now fi is replaced by a time-dependent which follows a multivariate random walk model, 
known also as local level model (Durbin and Koopman 31 ). 

Discount Weighted Regression (DWR), which originated in the path-breaking work of 
Brown 32 , is a method for forecasting autocorrelated time series. Considering univariate time 
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series Ameen and Harrison 33 developed further DWR for more complex time series. The 
reviews of Ameen 34 and Goodwin 35 suggest that DWR can model efficiently time series in a 
wide range of situations. Triantafyllopoulos and Pikoulas 36 developed a multivariate version 
of DWR and these authors focused on the estimation of the measurement covariance matrix. 
In this paper we consider the DWR method of Triantafyllopoulos 37 for multivariate local level 
models defined by 

yt = lit + e t and yn = Ht-i + Ut, (2) 

where ej ~ A/^(0,£) and u>t ~ Ap(0, J)tS). The scalar is specified with the aid of a 
discount factor 5 and the sequences {ef } and {u)t\ are mutually and individually uncorrelated, 
e.g. E(eieJ-) = E(uj k u' e ) = E(e r u:' s ) = 0, for all i + j, k ± I and for all r,s. Here E(-) 
denotes expectation and e'j denotes the row vector of €j . The model definition is complete by 
specifying a prior distribution p(fio\Y,), which is usually the p-dimensional normal distribution, 
e.g. jUo|S ~ N p (rriQ, PqT,), for some known prior mean vector mo and a positive scalar Pq > 0. 
It is further assumed that fiQ is uncorrelated of all w;. For some positive integer N > 0, let 
U = (yiiVZi ■ ■ ■ iVt) be the information set comprising data up to and including time t, for 
t = l,2,...,N. 

With the prior fio\T, ~ M p (mo, PoT>), the posterior density of /xj|E, y* is /^|S, y* ~ 
M p (m t , -PfS), where m t and P t are updated by 

Pt-x Smt-i + Pt-Wt j D 1 f o\ 

mt = mt - 1 + TTp7 l et= 5 + and Pt = JTpT^ (3) 

with et = yt — E(y t |y* _1 ) = y t — mt-i being the one-step forecast error vector at time t — 1. 
Define the residual error vector r t = E(et|y*) = yt — m t . For each time t the estimator St of 
S is achieved by least squares estimation as 



after observing that 



r t = y t - m t = y t - mt-i - ; , „ — = 



-Pt-iet i^-iet fe 



5 + P t _i <5 + Pt_i 5 + Pt-x 

Details of the derivations of mt, Pt and St appear in Triantafyllopoulos and Pikoulas 36 and 
Triantafyllopoulos 3 . 

From the above it follows that the one-step forecast density is 

t f (8 + P t )St 
yt + i\T, = b t ,y ~N p <m t , 



and the corresponding one-step forecast error density is 



et + ^ = S t ,y^ K ^l±p^l^ 



(5) 



where e t+ i = y t +i - E(yt + i|y 4 ) = y t+1 - m t . 

The adequacy of the model is evaluated via the mean of squared standard one-step forecast 
error vector (MSSE), the mean of absolute percentage one-step forecast error vector (MAPE) 
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and the mean of absolute one-step forecast error vector (MAE). These statistics are discussed 
in Chatfield 38 and for data y±, 1/2, • • • ,Vn they are defined by 

1 N 

MSSE = -J2 Mtf (eltf ■ ■ ■ (eU?]', 4 = 
t=i 

mae = ^Eo e «i i e2 *i ••• W> 

t=i 

where e* is the standard one-step forecast error, y t = [yu yit • • • Vpt]', &t = [eit e 2t • ■ • ept]' 
and {5~ 1 (5 + Pt~\)St-i\~ 1 ^ 2 denotes the inverse of the symmetric square root of the matrix 
5~ 1 (5 + Pt-i)St-i based on the spectral decomposition of symmetric matrices (Gupta and 
Nagar 39 ; pages 6-7). If the model fit is good the MSSE should be close to the vector [11 • • • 1]', 
while MAPE and MAE should be as small as possible in absolute value. Note that the MAPE, 
as a percentage statistic, makes sense only for a positive valued process yt, for all t. If this is 
not the case, then MAPE can not have a meaningful interpretation and it should be excluded 
from the statistical analysis (Chatfield 38 ). 



{ s } e > 



MAPE = — V 



4- 1 



\ e u\ \e-2t\ 



yu y2t 



M 

Vpt J 



3 The Bayesian Control Chart 
3.1 The Main Idea 

Bayes' factors have been extensively discussed in the statistics literature and recently they 
have been applied sequentially for time series, see e.g. West and Harrison 28 (Chapter 11). 
Salvador and Gargallo 40 propose a monitoring scheme, based on Bayes' factors, for multivari- 
ate time series, but this approach is not suitable for control charting, because it is applied in 
a model selection problem. In addition to this, most of the Bayesian time series monitoring 
(including the work of Salvador and Gargallo 40 ) relies upon simulated based methods and in 
particular Monte Carlo simulation. In this paper we favour non-iterative techniques, because 
they are faster, more flexible and easier to apply. 

Once we have the distribution ([5]) we can construct a target distribution for the dispersion 
of yt from the target mean and then compare these two distributions. It is well known (see e.g. 
Pan and Jarrett ) that the forecast errors and ej (i 7^ j) are approximately uncorrelated 
and the approximation is so good as St is closer to S. Suppose now that the target mean 
of {yt} is denoted by /i and the process dispersion covariance matrix is denoted by V. This 
notation is consistent with the Shewhart-Deming model as in equation ([1]), with V = T, so that 
E(yj) = fj, and Var (yt) = S, where Var (yt) denotes the covariance matrix of yt- Is is assumed 
that fj, is a generally unknown vector, but not stochastic. In our model of equation ([2]) we 
have K(yt\fit) = IH an d Var (7^ | /if) = S, but now fit is stochastic and it also changes with time 
according to the random walk model of ([2]). We postulate that, if the process is in control, 
the one step forecast mean of yt will be close to the target mean vector /i and the forecast 
covariance matrix of yt will be close to the target dispersion covariance matrix V. Thus we 
can define the target error distribution by e t ~ N p (0, V), where Et = yt — fJ> is the process error, 
also known in the process adjustment literature (del Castillo 3 ) as disturbance drift. Here we 
assume that V is positive definite matrix, although the proposed approach can be modified 
when V is positive semi-definite. According to the above postulate, if model (H|) describes 
well the in-control process, density ([5]) should be close to the above target distribution. In 
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order to find out "how close" it is, we form the Bayes' factor at time t: 

BF{t) -m- — m — ' t - 1 ' 2 '--- iV ' 

where f e (t) and f £ (t) denote the probability density functions of et and £t, respectively. 

For consistency in the above equation we need to make the convention y° = (the null or 
empty set). Since both densities f e (t) and f e (t) are normal we have 



-S(y, - mi-O'S^ibd - m,_i)/(M + 2F t _0} , (6) 

where det(-) denotes the determinant of a square matrix. The Bayes' factor BF(t) takes 
values from to +00. We will say that the process yt is in control at time t, if f e (t) = f e (t), 
or if BF(t) = 1; otherwise the process will be out of control, at this time point. An out 
of control signal might be caused because of a mean shift (e.g. when E(y<|y* _1 ) = mt—i is 
significantly different than //) or because of a dispersion shift (e.g. Var(y t |S = St_i,y ) = 
(5 + Pt-i)St-i/5 is significantly different than V). 

3.2 The Modified EWMA Control Chart for Correlated Data 

A control chart for the Bayes' factor BF(t) can conclude whether BF(t) is close to 1 and thus 
whether the process is in control or not. Since BF(t) is positive valued, it is more convenient 
to work with the logarithm of the Bayes' factor 

LBF(t) = log BF(t) =plog(5/2 + {logdet(F)}/2-p{log(<5 + J P t _i)}/2- 
{logdet(5 t _ 1 )}/2 + (yt - ri'V- l (y t - fi)/2 

-S(y t - mt-iYS-^yt - m t _i)/(2<5 + 2P t _i) (7) 

and so we can construct an appropriate univariate control chart for LBF(t). In order to 
propose such a chart we need to deal with two issues: (a) the values of LBF(t) will be 
serially correlated and (b) the distribution of LBF(t) might not be normal. 

Considering (a), in our development it is clear that, from the definition of the BF(t), either 
the original data yt are i.i.d. or auto-correlated, the resulting data BF(t) (or LBF(t)) will be 
correlated and hence, if the Shewhart or any other control chart is to be used successfully, they 
should be modified appropriately to accommodate for correlated observations. Many authors 
have demonstrated that the Shewhart control charts need to be modified in order to cater 
for serially correlated observations (Vasilopoulos and Stamboulis 41 ; Schmid 42 ). Similarly, the 
EWMA needs also to be modified and the resulting modified EWMA control chart has been 
discussed in many articles including Schmid and VanBrackle and Reynolds 44 . According 
to Harris and Ross 43 ignoring serial correlation has a stronger effect in EWMA than in the 
Shewhart control chart, but as we will see later the EWMA control chart is preferable to 
Shewhart, because it is more robust to the assumption of normality. One could also consider 
the modified CUSUM chart for correlated observations, but we will not further discuss this 
in the present paper. 

Proceeding with (b) one needs to check the assumption of normality, before applying a 
modified EWMA (or Shewhart or CUSUM) control chart. Borror et al. studied the ARL 
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performance of the EWMA and they suggested that the EWMA with a smoothing parameter 
equal to 0.05 is very effective, even in the presence of non-normality of the observations. 
This result agrees with Montogomery 1 who states for the EWMA "It is almost a perfectly 
non-parametric (distribution free) procedure". Maravelakis et al. 47 study the robustness to 
normality of the EWMA by tabulating characteristics of the run length distributions (e.g. 
ARL) for observations generated by several gamma distributions. These results conclude 
that, for relatively low values of the damping parameter of the EWMA and for shifts in the 
mean the EWMA control chart can be used, even in the absence of normality. Moreover, 
if the process is in-control following a symmetrical, but not normal, distribution, then the 
EWMA can be applied successfully. To the following we look at the empirical distribution of 
LBF{t) when the process is in control and when it is out of control. 

We generate 1000 vectors from a bivariate normal distribution A/2(/f, V) with 








and V 



1 2 

2 5 



and we generate 1000 vectors for three out of control scenarios. In scenario 1 we simulate data 
from A/2 (/id) V) (deviations from the mean /i); in scenario 2 we simulate data from A/"2(/i, Vd) 
(deviations from the covariance matrix V); in scenario 3 we simulate data from A/2 (/id, Vd) 
(deviations from both /t and V), where 



/id 



0.5 




and Vh 



1 

2.5 



2.5 



Figure [T] shows the histograms of the LBF(t) for the above four scenarios (one in control and 
three out of control scenarios). From this figure we observe that, although the distribution 
of the LBF(t) for the in-control process (panel (a) in Figure [1]) is not-normal, it is roughly 
symmetric. The distributions of the LBF{t) for the out of control processes appear to be 
slightly skewed, but the histograms are not conclusive. The important point is the non- 
normality of the LBF{t) and the symmetry of the distribution of the in-control process. 
This enables us to make use of the modified EWMA control chart, but we note that the 
modified CUSUM control chart can also be used. A more formal confirmation of the non- 
normality of the distribution of LBF(t) can be carried out by the using standard tests of 
normality, however, here the histograms are deemed sufficient to declare the non-normality 
of the distribution of LBF(t). 

We use a two phase control scheme; in Phase I the mean /t and the covariance matrix £ 
are estimated and adjustments are applied if necessary, while in Phase II the EWMA control 
chart is applied to detect any changes in the mean of LBF{t). Thus we propose the algorithm: 

Algorithm 1. There are two phases: 

Phase I: We fit the DWR model (0|) for a set of historical data t = 1, 2, . . . , N* , with N* < 
N . We check the performance and adequacy of the model via the MSSE, MAPE and 
MAE over all t = 1,2, .. . , N* and we possible apply adjustments to the DWR model, 
(e.g. adjustments in the mean level) so that we obtain optimal values m op t = tun* , 
Sopt — Sn* , <5 = 5 apt ensuring that in Phase I the model matches the in-control process. 
The modified EWMA control chart is applied so that control limits are adequately defined 
according to pre-specified ARL curves. For this to be designed, a state-space model for 
the process LBF{t) needs to be identified and here simple AR and ARMA modelling will 
be generally acceptable. 
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Figure 1: Histograms of the log Bayes' factor LBF(t) for an in-control process (panel (a)) 
and out-of-control processes (panels (b)-(d)). The out of control scenarios considered are 
deviations from the mean vector (panel (b)), deviations from the covariance matrix (panel 
(c)) and deviations from both the mean vector and the covariance matrix (panel (d)). 

Phase II: We fit the DWR model with the model components from Phase I (e.g. 5 = 5 op t, 
m t = m opt, S = S op t and we apply a modified EWMA control chart at observations 
LBF(t) with the control limits identified at Phase I, for t = N* + 1, N* + 2, . . . , N. 

In order to apply the modified EWMA control chart we first calculate the series zt with 
observations xt = LBF(t) as 

z t = Xx t + (1- X)z t . 1 , 0<A<1. (8) 

The parameter A is the EWMA smoothing parameter and as it is mentioned above, for 
A = 0.05 or A = 0.1 the control chart is robust to normality. Then, the control limits of the 
modified EWMA control chart are 

fi z ±ca z , (9) 

where fj, z = E(z t ), a\ = lim^oo Var(z t ) (asymptotic variance of zt) and c > is determined 
according to the required ARL. For AR(1) dependence xt = 4>xt-i + v% and for large t, the 
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asymptotic variance a 2 is 




a 2 A{l + 0(l-A)} 



(l-02)(2-A){l-0(l-A)} 



where v t ~ A/"(0, a 2 ) and a 2 , are assumed known. In practice these parameters are estimated 
at Phase I. According to Schmid 43 the asymptotic variance a 2 performs better than the exact 
variance of zt, which is given in Schmid 43 and which produces time-dependent control limits. 
Most of the literature on this topic focuses on deriving the variance a 2 assuming simple time 
series models for xt, e.g. as in the above AR(1) or as in the ARMA(1,1) model considered in 
VanBrackle and Reynolds 44 . 

Algorithm [1] can be simplified, if at Phase I, the quantities Pt and St converge to stable 
values and these values are determined in Phase I for both phases. This brings up a well 
known problem, which has received considerable attention in the time series literature (see 
e.g. Durbin and Koopman 31 ). However, for the DWR and similar multivariate models limiting 
results for Pt and St have not been yet established. The next theorem (which proof is in the 
appendix) states that Pt and St converge to stable limiting values. 

Theorem 1. In the DWR model (0|) the estimator St of the measurement covariance matrix 
£ converges in probability to £ and the non-stochastic scalar parameter Pt converges to the 
limit P = (V<5 2 + 4 - 8)/2, i.e. S t £ and P t — ► P. 

From Theorem[T]the estimator St is consistent and from the proof of this theorem (given in 
the appendix), St is also unbiased estimator. Theorem Q] suggests that Pt-\ in the calculation 
of LBF(t) of equation ([7]) can be replaced by its limit P. From equation ([3]) and Theorem Q3 
the forecast of yt, mt-i can be approximated by 



where Pt-i of equation ([7|) is replaced by P. Figure [2] shows how fast {Pt} converges to its 
limit P, for a prior Pq = 1/1000 and three values of 5. This figure points out that Pt is 
bounded above by 1, but for 5 = 0.2, this bound is only achieved after t > 13 (solid line in 
Figure [2]), while for 5 = 0.9, this bound is achieved for any t > 1 (dotted line in Figure [2]). 
This gives an empirical indication of the speed of convergence of {Pt}, for several values of 5. 

The limit P is known before the algorithm starts (e.g. P depends only on 5) and, given 
enough data in Phase I, the limit £ can be approximated by E ~ SV* , in the end of Phase 
I. This can have an additional benefit on computational savings, but more importantly it 
gives a theoretical justification that the DWR produces a good copy of the process {yt} and 
therefore this model is appropriate for the monitoring part at Phase II of Algorithm [TJ For 
example, if Pt and St were not converging to stable values, no matter how many data we 
collected at Phase I, the covariance matrix of yt and thus its uncertainty would change over 
time resulting in an unstable time series model. False alarms are probable in the framework 
of such unstable models, which should be avoided. 

In the design and application of the control chart it is important to suggest values of mo, 
Pq, 5 and So an d to study their sensitivity and influence to the performance of the proposed 
control chart. Since these suggestions are related to forecasting as in equation ©, results 
on the sensitivity of such prior parameters follow from Triantafyllopoulos and Pikoulas 36 and 
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Figure 2: Rate of convergence for the sequence {Pt} of Theorem [TJ the solid line plots {Pt} 
for 5 = 0.2, the dashed line plots {Pt} for 5 = 0.5, the dotted line plots {Pt} for 5 = 0.9 and 
the dashed/dotted line is the critical bound of 1. 

Triantafyllopoulos 37 . It is worthwhile noting that, given enough data in Phase I, the values of 
mo, Po and So are not critical to the forecast performance, as in time series modelling prior 
information is deflated over time. This is indicated in Theorem Q] from the fact that P does 
not depend on Po- The value of 5 can be critical in forecasting and a general recommendation 
is that several values of S (in the range of (0, 1)) are applied in Phase I and according to the 
forecast performance (see Section [2]) a value of 5 is decided. One should note that high values 
of 5 (e.g. 5 = 0.9) yield smooth forecasts with low forecast variances, but these forecasts are 
sometimes unable to forecast abrupt changes in the data; low values of 5 (e.g. 5 = 0.1) yield 
more precise forecasts in the presence of "wild data" , but these forecasts come with increased 
forecast variances. 

Our proposal for the modified EWMA control chart for the LBF(t) process is motivated 
from the fact that the observations LBF(t) possess autocorrelation and non-normality. The 
approach is model-based, and so a comparison with traditionally used multivariate control 
charts, such as the Hotelling's T 2 and the M-EWMA (which are both data-based control 
charts), is difficult and in many occasions it can not give justice. Within the model-based 
control charting methods, it appears that our approach can be compared with the residual 
chart (Pan and Jarrett 24 ), but again the comparisons need to make sure that model uncer- 
tainty (whether for example the DWR is a good model or an alternative time series model 
performs better) should be ideally removed before any comparison is attempted. For example 
a miss-specification of a time series model might result to a false result in the comparison 
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of the competing control charts. From our experience the DWR works generally well (since 
it is a generalization of the Shewhart-Deming model), but this might not be the case for 
every multivariate process. We believe that such a comparison should deserve the length and 
the detail of a whole paper and thus here we do not pursue this project. Next we give two 
examples illustrating the design and application of the proposed control chart. 

To the above we have assumed that given a process {yt} the interest is in building a 
control chart for monitoring simultaneously the process mean and the dispersion covariance 
matrix. However, in some cases the interest is placed on monitoring the dispersion covariance 
matrix only. In this case we can modify the control scheme by considering a modified EWMA 
control chart of the log-Bayes' factors of the first order difference process zt = yt — Vt-u 
which from equation ([2]) has zero mean. Control charts based on {z{\ will be more robust as 
compared to those for {yt}, since the uncertainty of monitoring the process mean of {yt} has 
been removed. 



4 London Metal Exchange Data 

London metal exchange (LME) is the world's premier non-ferrous metals market trading 
currently aluminium, copper, lead and zinc, among other non-ferrous metals. Information on 



the LME and its functions can be found in its web site: http : / / www . lme . co . uk| The review 
of Watkins and McAleer 48 explores the recently growing literature on the LME market and 
Triantafyllopoulos 37 discusses the correlation of spot and future contract prices of aluminium 
based on the DWR model of Section [2j In this paper we discuss data of spot prices for the 
four metals aluminium (variable {yit}), copper (variable {y2t}), lead (variable {?/3t}) and zinc 
(variable {yu})- 

The data are collected from January 2005 until October 2005 for every trading day ex- 
cluding weekends and bank holidays; Figure [3] plots the data. We form the observation vector 
Vt — [yit U2t U3t Vit]' and we are interested in knowing whether volatility is apparent, for 
t = 151 until t = 220. In other words we want to know whether from t to t + 1, the variability 
of the observations yt and yt+i has changed. This is a major concern to econometricians, 
because if there is evidence for volatility, this means there is uncertainty in investments and 
ideally the volatility should be understood and explained. In order to answer this important 
question we form the first order difference of the series {yt}, defined by xt = yt — yt—i, for 
t > 1 (Figure HJ) . Adopting the usual forecasting strategy of commodity forecasting, given 
data up to time t — 1, the forecast mean of yt at time t is just the value of yt-i and so 
we can write E(yt|y* -1 ) = yt-\- We note that the true mean of xt may not be zero (unless 
in model ([2]) it is = fj, + u>t), but it is true that conditionally on y t_1 or yt-i we have 
^( x t\y t ^ 1 ) = lE(y<) — yt-i = 0, since K(yt\y t ~ 1 ) = yt-i- From Figured] we observe that the 
series {xt} fluctuates around zero and volatility can be detected as significant deviations from 
the zero target; such deviations can be detected with the aid of a control chart of Section [3j 

First we need to make sure that the DWR model fits the differenced series {xt} well. We 
take t = 1 — 150 as Phase I, in which the adequacy of the DWR model is evaluated. The 
performance statistics of Section [2] are: MSSE = [0.993 1.486 0.866 1.323]' and MAE = 
[18.932 45.187 14.569 19.082]', suggesting an acceptable fit. Of course the MAPE is not 
available, since {xt} is not a positive valued process (Section [2]). 

We have designed a modified EWMA control chart for the LBF{t) of the process {x{} 
according to the discussion of Section Figure [5] shows four control charts corresponding to 
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London metal exchange data 




Time 

Figure 3: LME data y t = [y lt y2t Vzt Vit]' , consisting of aluminium ({y u }), copper ({y 2 <}), 
lead ({y3t} and zinc ({yu}) spot prices (in US dollars per tonne of each metal). 

four values of the EWMA smoothing parameter A. Typically the control chart is robust to 
normality for small values of A, but for these values the control chart is only detecting very 
small drifts in the mean this might not be desirable. As A increases the modified EWMA 
control chart is losing its robustness over normality, but for symmetric process distributions, 
such as the empirical distribution of the LBF{t) shown in Figure [U the EWMA control 
chart might still be used for A = 0.5. The correlation of the LBF{t) is accounted by the 
autoregressive model of Section [3] and an analysis involving the data at Phase I shows that an 
the autoregressive parameter (ft = 0.1 is adequate to capture the autocorrelation of LBF{t). 
According to Tables for the ARL of the modified EWMA control chart (see e.g. Shiau and 
Hsu ) we choose the value of c in equation (|9|) so that ARL = 370.4, e.g. for A = 0.05 and 
eft = 0.1 we have c = 2.469. The remainder of the control limits are calculated as in equation 

©■ 

Figure [5] shows that the process in Phase II appears to be in control, for A = 0.05 and 
A = 0.1, while for A = 0.2 and A = 0.5 the control chart returns an out of control point 
at t = 172 (with values zm = —1.852 and zm = —2.999, respectively). The mean of the 
EWMA zt is slightly lower than zero, which indicates that, for the entire process {xt}, there 
will be some deviation of the predictive density f e from the target density f e . It is up to the 
modeller to decide whether such deviation from the target distribution is worth of declaring 
the process out of control. In search of a more automatic approach, one can lift up the 
whole control chart so that in Phase I the mean of Zt is exactly zero. This can be performed 
automatically, in the end of Phase I, and this will declare the process in control in Phase II, 
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London metal exchange differenced process 




Time 

Figure 4: LME differenced process xt = [x\t Xit x^ t x^t]', consisting of aluminium ({xit}), 
copper ({x2t}), lead ({x^t}) and zinc ({xm}). The horizontal lines, placed at zero, indicate 
no volatility. 

for A = 0.05,0.1, while for A = 0.2, A = 0.5 there is an out of control point at t = 172. In 
Figure [5] the value of A = 0.5 is rather high to ensuring correct control limits of the modified 
EWMA chart (see the relevant discussion in page [7]); here the chart with A = 0.5 is mainly 
shown for comparison purposes with the charts with lower values of A, but in practice we 
suggest that A does not exceed 0.2, unless there is strong evidence to support the assumption 
of normality for the distribution of LBF(t). It is worth pointing out that the concentration of 
consecutive EWMA values under the mean in Phase II is causing warning, which is apparent 
in all charts. The phenomenon is more apparent in the charts for A = 0.05 and A = 0.1 and it 
can suggest the out of control state of the process at t = 172, which is apparent in the charts 
with A = 0.2 and A = 0.5. The interpretation of the out of control signal at t = 172 can not 
be done just by looking at Figure H] and more dedicated methods of out of control variable 
identification need to be employed, see e.g. Bersimis et al. 14 . 

5 Production Time Series Data 

In an experiment of production of a plastic mould the quality is centered on the control of 
temperature and its variation. For this purpose five measurements of the temperature of 
the mould have been taken, for 276 time points. The experiment is fully described in Pan 
and Jarrett 24 and these authors show that this 5-dimensional production process {yt} is both 
autocorrelated and serially correlated including both vector autoregressive and moving aver- 
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(a) lambda = 0.05 (b) lambda = 0-1 
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Figure 5: Modified EWMA control chart for the log Bayes' factor of the LME differenced 
process. Plots (a)-(d) show the modified control chart for different values of the smoothing 
parameter A. In each plot of the panel, the solid horizontal line indicates the mean of the 
EWMA and the dotted horizontal lines indicate the control limits; the vertical line separates 
Phase I (for t = 1 - 150) and Phase II (for t = 151 - 210). 



age terms. These authors use a vector state space charting approach based on the Hotelling 
control chart resulting on 12 out of control signals at Phase II (time points from t = 181 to 
t = 220) and hence concluding that the process falls badly out of control at Phase II. 

We have used the data at Phase I (time points t = 1 — 180) in order to estimate the 
target mean vector /i = [208.245 153.638 53.063 —22.742 16.126]' (as the average of each 
t = 1 — 180) and the dispersion covariance matrix 



V 



0.168 -0.001 

-0.001 0.023 

0.633 -0.017 

-0.438 0.006 

0.015 -0.002 



0.633 -0.438 0.015 

-0.017 0.006 -0.002 

25.621 -15.658 0.453 

-15.658 14.181 -0.596 

0.453 -0.596 0.951 



(as the sample covariance matrix of each y t : t = 1 : 180), where yt = [yit, V2t,V3t, VAt^Vbt]' ■ The 
DWR fits well with MSSE = [0.855 0.950 0.992 1.161 0.996]', which is close to [1 1 1 1]. The 
other two performance statistics are MAE = [1.378 0.899 4.450 3.316 0.945]' and MAPE = 
[0.007 0.006 0.089 - 0.059]', where for {y 4t } the "-" indicates that the MAPE is not available, 
since this variable is not positive valued (see the relevant discussion for MAPE in Section [2]). 
The above performance statistics suggest that the model fit is good and therefore we can 
proceed with control charting at Phase II {t = 181 — 279). 

The first thing to do is to find a suitable AR(1) model for the process LBF{t). A suitable 
model is the AR(1): LBF{t) = -4.624 + 0.062LBF(t - 1) + v t . According to the discussion 
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(a) lambda = 0.05 (b) lambda = 0.1 

Figure 6: Modified EWMA control chart for the log Bayes' factor of the Production process. 
Plots (a)-(b) show two charts for values of the smoothing parameter A = 0.05 and A = 0.1. 
For both plots, the solid horizontal line indicates the target mean and the dotted horizontal 
lines indicate the control limits; the solid vertical line separates Phase I (for t = 1 — 180) and 
Phase II (for t = 181 - 276). 



above, we remove the intercept —4.624 so that we can obtain a in-control process in Phase I. 
Thus we design the modified EWMA control chart for LBF{t) + 4.624. Again we use tables 
for the modified EWMA control chart and for A = 0.05 the resulting control chart is given 
in Figure El This figure agrees with the residual chart of Pan and Jarrett 24 , that finds the 
process in Phase II out of control for most of the data points. In Phase I chart of panel (b) 
of Figure [6] gives one out of control point, which is in agreement with Pan and Jarrett 24 , but 
in panel (a) of Figure [6] the control chart detects more out of control points in Phase I. The 
EWMA control chart is robust to non-normality for the low values of A = 0.05 and A = 0.1, 
but for A = 0.05 the chart is more sensitive to small shifts in the mean of LBF(t), resulting to 
the detection of out of control points in Phase I. Any out of control points in Phase I should 
be immediately investigated and usual SPC procedures of removing influence of these points 
in the calculation of the control limits should be applied (Montgomery 1 ). 



6 Conclusions 

This paper develops a new multivariate control chart based on Bayes' factors. This control 
chart is specifically aimed at multivariate autocorrelated and serially correlated processes. 
The general idea is to form a target distribution, to construct a predictive density with 
good forecast ability and then to apply a univariate control chart for the logarithm of the 
Bayes' factor of the predictive error density against the target error density. Although in this 
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paper, for simplicity, we have considered normal distributions for the target and the predictive 
densities, in general application the proposed control charts can be applied considering other 
densities too as long as they are available in analytic form. 

We have restricted our discussion to the modified EWMA control chart, but other control 
charts such as the modified CUSUM and non-parametric control charts can be applied. A 
major advantage of our approach as compared to other multivariate control charts is that once 
we have obtained the log Bayes' factors we can apply any appropriate univariate control chart. 
A difficulty appears to be that the resulting Bayes' factors process is both autocorrelated and 
non-normal, but we believe the design of the proposed chart is a challenge that can attract 
and motivate further research in this so important area of statistical process control. 
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Proof of Theorem^ First we prove St — — ► E. It suffices to prove that St is unbiased estimator 
and that its covariance matrix converges to zero. From equations (j4]) and ([5]) we obtain 



8(5 + P,_!)S 1 



<5 + P;_i (S + P^S t 



(iE) 



and so St is unbiased for E. For the convergence, let vech(-) denote the column stacking 
operator of a lower portion of a covariance matrix and let || • || denote a matrix norm defined 
in a suitable linear space. From equation ([5]) we have 



1 1 

Var{vech(5 t )} = ^^ 



5 + Pi- 



Var{vech(eje£)}. 



(A-1) 



From equation ([5]) ej follows ap-variate normal distribution and so by writing ej = [en ea ■ ■ ■ e 
we have that Cov(ejj, e^) = E^^-e^) are bounded, since these expectations are expressed as 
moments of the multivariate normal distribution (Triantafyllopoulos 50 ). Hence Var{vech(eje-)} 
has finite elements and so we can write || Var{vech(eje£)} ||< M, for some M > 0. For any 
e > define t = [eM] (the integral part of eM). From Pj_i > we have that 6/(6+Pi-i) < 1, 
for all i = 1,2, ...t. Then 



'■Pi ; 



||Var{vech(S t )}|| 



1 



M 

5 F 



E 

i=i 
t 

E 



8 + Pi-i 



Var{vech(eje^)} 



<5 + P, 
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tM M 

* ^ = T <e ' 



for any t > to- This shows that limj^oo Var{vech(S'f)} = and so St — > E. 
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Proceeding now with {Pt} we show that {Pt} is a Cauchy sequence in the real line and 
hence lim^oo Pt = P exists. To prove that {Pt} is a Cauchy sequence, it suffices to prove 
that linn^oo \P t — Pt-i\ = 0, where | • | denotes absolute value. First we show that exists 
positive integer to such that for all t > to it is Pt < 1. The proof of this is by contradiction. 
Suppose that for all to exists t > to such that Pt > 1. Without loss in generality take to = t* 
and P t * = 1. Then we see that P t * +1 = 1/(5 + P t .) = 1/(5 + 1) < 1, P t * +2 = 1/(5 + P t *+i) = 
(5 + l)/(5 2 + 5 + 1) < 1 and likewise P t *+k < 1, f° r a U k > 1. So we can pick t = t* + 1 so 
that we can not find any t > to with Pt > 1, which contradicts the hypothesis. Thus exists 
to > so that for all t > it is Pt < 1. This in turn implies that 

5 + P t -i>l, Vt>t . (A-2) 
From the definition of P( of equation ([3]), we obtain 



1 1 Pt-x-Pt-i (-l) t - l (P 1 -P [ 



0J 



P * P< 1 * + P-i * + P-2 (5 + Pt_ 2 )(5 + P t _ 2 ) "' 1^1(5 + P t -i)(5 + P M ) 

Now pick t as in (|A-2j) and define M = min{5 + P t _i, (5 + P*_2) 2 , ...,(£ + Pi +i) 2 } so that 
M > 1. Then 



IP_ P , ji-jft-gj ii-jpq-Pq 2 ! 



n:io(^ + p) 2 ritr^ + p-i)(* + p-*-i) nr=o(«^ + p) 2 m<- 

since lim^oo M (_to_1 = +00. This proves that lim^oo — Pt— 1 1 = and so {Pt} is a Cauchy 
sequence. Thus lim^oo Pt = P exists and from equation ([3|) we have P = l/(5 + P), for which 
we derive P = (VS 2 + 4 — £)/2, after rejecting the negative root P = (— \J 5 2 + 4 — <5)/2. □ 
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